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A THEORETICAL EXAMINATION OF DIFFUSIVE 
MOLECULAR DYNAMICS 

G. SIMPSON, M. LUSKIN, AND D.J. SROLOVITZ 


Abstract. Diffusive molecular dynamics is a novel model for materials 
with atomistic resolution that can reach diffusive time scales. The main 
ideas of diffusive molecular dynamics are to first minimize an approxi¬ 
mate variational Gaussian free energy of the system with respect to the 
mean atomic coordinates (averaging over many vibrational periods), and 
to then to perform a diffusive step where atoms and vacancies (or two 
species in a binary alloy) flow on a diffusive time scale via a master equa¬ 
tion. We present a mathematical framework for studying this algorithm 
based upon relative entropy, or Kullback-Leibler divergence. This adds 
flexibility in how the algorithm is implemented and interpreted. We 
then compare our formulation, relying on relative entropy and absolute 
continuity of measures, to existing formulations. The main difference 
amongst the equations appears in a model for vacancy diffusion, where 
additional entropic terms appear in our development. 


1. Introduction 

One of the outstanding challenges in atomistic simulation of condensed 
systems, such as solids, liquids, and glasses, is access to experimentally 
meaningful length and time scales. The spatial scales amenable to direct 
molecular dynamics (MD) simulations have grown over time; MD simu¬ 
lations on current large scale parallel computational facilities now reach 
over 10^^ atoms or approximately 1 //m^, [7]. Significantly larger length 
scales can be achieved by application of multiscale modeling techniques 
that combine atomistic simulations with continuum methods (see, for ex¬ 
ample, [Il[2l[T6l[22l[26l[28l[33]). With regard to time, MD simulations have 
fundamental time scales associated with atomic vibration periods (~ 10“^^ 
s), but typical MD time steps are two orders of magnitude smaller. The 
longest times that have been achieved in large scale MD simulations on 
special purpose hardware is ~ 10“^ s, |30] . More typically MD simulations 
access times of less than than 10“® s. Hence, reaching laboratory time scales 
remains amongst the most important challenges in the application of MD 
today. 

Many approaches have been developed to address the molecular dynamics 
time-scale challenge. Of these, methods based upon transition state theory 

Date: June 9, 2015. 

GS and ML were supported by US Department of Energy Award DE-SG0012733. 

1 



2 


G. SIMPSON, M. LUSKIN, AND D.J. SROLOVITZ 


have been widely applied and their continued development remains an ac¬ 
tive area of research. Transition state theory-based methods rely on rare 
event ideas in which the system explores a particular energy basin for a 
long period of time and makes infrequent transitions from one basin to an¬ 
other, [9l|T3l[23l3Tl[32] . Access to long time scales is provided by replacing 
the true dynamics within the basin with a stochastic approximation and 
the prediction of the times between the transitions to other basins. Such 
methods rely on the transitions being sufficiently rare; i.e., the computa¬ 
tional time required to characterize the dynamics within the basin is much 
shorter than the time between inter-basin transits. The main difficulty with 
this approach arises when the transitions between basins are not rare; for 
example, at high temperature or in situations where there are low energy 
barriers between basins. Another challenge arises in very large systems as 
the time between rare events occurring somewhere in the system typically 
scales inversely with system size. Nonetheless, progress has recently been 
made to extend these methods to large systems, [TTlfT^IM] . 

An important distinction between classes of “events” in many materials 
science applications is between those that are displacive and those that are 
diffusive. Displacive events are often collective and involve relative atomic 
displacements that are small compared with a typical interatomic separa¬ 
tion (~ 2 A). Diffusive events, on the other hand, often involve a series of 
atom or vacancy hops amongst atomic sites. In a solid, the time between 
hops is commonly many orders of magnitude larger than the atomic vibra¬ 
tion period. This type of scale separation is necessary for the application 
of traditional transition state theory-based approaches. The motion of de¬ 
fects in materials and many types of phase transformations occur through 
a combination of displacive and diffusive events. 


1.1. Diffusive Molecular Dynamics. One method that explicitly takes 
advantage of this separation in scale between diffusive and displacive events 
in atomistic simulations is the so-called Diffusive Molecular Dynamics (DMD) 
method |15 t l24|. The idea is to first minimize an approximate free energy 
(including atomic bonding and atomic vibrational degrees of freedom) of the 
system with respect to the mean atomic coordinates (averaging over many 
vibrational periods), and then to perform a diffusive step where atoms and 
vacancies ffow on a diffusive time scale. This free energy is typically de¬ 
scribed using the variational Gaussian (VG) method [141118] . The time step 
for such simulations is the relatively long (compared with vibrational) diffu- 
sional time scale. In this approach, the search for transition state barriers is 
replaced by the introduction of an empirical diffusion coefficient or mobility. 
While this coefficient may be determined directly from atomistic simulation 
for every local transition, it is commonly viewed as a constant across the en¬ 
tire simulation; this constant can also be determined from purely atomistic 
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calculations. This approach is physically sensible when most of the impor¬ 
tant diffusive events are of the same type, as in an exchange of an atom for 
a vacancy. 

A key element of DMD is that, in contrast to traditional MD, each atomic 
site has associated with it a continuous probability of occupancy by an atom. 
See, for example, m- We denote this probability Cj at atomic site i. In the 
case of a single species (elemental) materials, Ci close to zero corresponds to 
site i being a vacancy with high probability, while Cj close to one corresponds 
to site i being occupied by an atom with high probability. Alloys can also 
be described in this manner. For example, q close to zero could correspond 
to species A while c* close to one could correspond to species B. This can 
be extended to multiple species, along with vacancies, by the introductions 
of additional degrees of freedom. Such an extension of the state space, to 
allow for longer time scale evolution, has also been explored in |29j . 


1.2. Relative Entropy. One of the tools we make use of in this analysis 
is relative entropy, or the Kullback-Leibler divergence (KL), [6]. KL is one 
of the many ways that the distance between two probability measures, such 
as ensembles, can be computed. It is broadly used in information theory, 
uncertainty quantification, statistical inverse problems, molecular dynamics, 
and other applications; see, for instance, [SlIinillZlEniEIlES] and references 
therein. 

Given two probability measures, u and /i on a common state space, the 
relative entropy distance between them is given by 


( 1 . 1 ) 


n{n\\u) 





^ -C u, 

otherwise. 


Here, p -C u if for any measurable set A such that u(A) = 0, we have 
/i(A) = 0 too. In this case, p is said to be absolutely continuous with 
respect to v. 

Relative entropy is a natural tool for this work as it can be directly con¬ 
nected to the Helmholtz free energy of a system. Indeed, the statement that 
7^ > 0 is equivalent to the Gibbs-Bogoliubov inequality [25]. To motivate 
this, consider the canonical ensembles associated with potential V{x) and 
an approximate potential V{x)-. 

(1.2) u{dx) = i){dx) = dx, 

where Z is the partition function, (3 = (/cbT)”^, /cb is the Boltzmann con¬ 
stant and T is the absolute temperature. Assuming that D <C 

(1.3) 7^(^>||u) = /3E^[H(x) - V{x)] -logZ + log Z. 

Dividing through by /3, and using Tl> 0, 

(1.4) -/3“^logZ</3E'^[AF]-/3-MogZ, 
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which is the Gibbs-Bogoliubov inequality. Mathematical treatments of TZ 
and its relationship to other popular measures of distance between proba¬ 
bility measures can be found in [6l|8]. 

1.3. Outline. In this work, we formulate DMD using relative entropy, to¬ 
wards the goal of constructing a rigorous mathematical foundation for the 
algorithm. We also examine under what assumptions the necessary ab¬ 
solute continuity will hold. This framework leads to naturally well posed 
variational problems that are inherently bounded from below by virtue of 
TZ being non-negative. 

Our paper is outlined as follows. In Section [2] we give our formulation of 
DMD using relative entropy. Then, in Section [3l we compare our expressions 
with those of the existing formulations of DMD. Several additional explicit 
calculations are given in the Appendix. 

2. A Relative Entropy Formulation of DMD 

To begin our description of DMD, we work in an extended state space 
that includes both the positions of the atom sites, x* € and the pres¬ 
ence/absence of atoms, Oj G {0,1}. For binary alloys, a* = 1 could corre¬ 
spond to species A and a* = 0 could correspond to species B. In either case, 
the total number of sites N is fixed and i = 1,... N. In what follows, we 
give separate developments of the vacancy and the binary alloy problems. 

2.1. DMD Potentials. We assume the existence of a potential V, that 
describes atomic bonding, as in traditional MD. Here, we focus on the class 
of pair potentials, 

(2.1) R(x) = -Xjl). 

i<j 

The associated DMD potential is 

(2.2) I/DMD(x,a) = ^ajaj(/>(|xj - Xjl). 

i<j 

In an elemental material, the potentials V and Vdmd are non-zero only when 
both atom i and atom j are present (i.e., the interaction between an atom 
and a vacancy is zero). In the case of a binary alloy (ignoring vacancies), 
we would have 

I/4R(x,a) = ^aiaj(pAA(lxi -Xj|) 

i<j 

(2.3) + + (1 “ ai)aj](pAB{\^i - Xj|) 

i<j 

+ ^(1 - ai){l - aj)(()ss(|xj - Xjl), 
i<j 

where (pTiT^ is the pair potential between a pair of atoms of types Ti and 
Tj G {A,B}. More sophisticated potentials, such as the embedded atom 
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method (EAM), can be used, but for the purpose of this work, it will suffice 
to consider the pair potential. 


2.2. Free Energy and Relative Entropy for the Vacancy Problem. 


2.2.1. Canonical Ensemble. Letting H be a bounded, open subset of we 
now formulate the canonical ensemble for a “cVT” (fixed total composi¬ 
tion, number of atomic sites and temperature) system, by introducing the 
distribution 

(2.4) v{dyi, a) = Z~^ exp {-/3(Vdmd - p • a)} dx, 

with /r = (pi,... jPat) is a set of chemical potentials. Chemical potential 
//j is associated with site i; each site has its own reservoir. While this is in 
the form of a generalized grand canonical ensemble, we view the chemical 
potentials as Lagrange multipliers enforcing the constraints 

(2.5) EZ[ai] = Ci, i = 

It is in this way that our v is akin to a canonical ensemble, with the mean 
occupancy values, the c, fixed. 

The partition function is 

(2.6) ^ / exp{-/3(VDMD(x,a) - p • a)}dx, 

and the ensemble averages require integration in space over D along with 
summation over all possible configurations of a = (oi,..., av). We thus 
assume: 


Assumption 1. For potential V, bounded set D, and occupancies c with 
Ci G (0,1), the chemical potentials pLi are finite and the partition function is 
finite and positive. 


We describe the distribution (j2.4l) as the “true” DMD distribution. The 
associated free energy is 

(2.7) E =-/?"Mog Z + p • c. 


In DMD, the concentrations are allowed to evolve under 


( 2 . 8 ) 


Ci = 


E 

j&N(i) 




\dcj 


dF 

da 


where kij is a mobility term, and N{i) includes neighbors of site i; both 
are, a priori, unconstrained. We note that (12.8j) is not the evolution studied 
by Sarkar [23], due to the challenge of computing the partition function; 
rather, they employ an approximate free energy. Nevertheless, we contend 
that, were computational complexity not an obstacle, (12.811 is the desired 
form of dynamics. 
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As the Cj evolve, the equilibrium positions of the atomic sites will also 
change, and can be obtained by computing 

(2.9) E"[x,]. 


2.2.2. Approximate Potential. Due to the computational challenge in eval¬ 
uating the partition function, (j2.6p . an approximate potential is introduced, 
V, with a corresponding probability distribution, u. This potential is then 
tuned to provide a “best” approximation of v by u. Here, we consider the 
following quadratic (harmonic) potential in the case of an elemental material 
(including vacancies) 

(2.10) HDMD(x,a;k,X) = 

i 

X and k are parameters approximating the mean atomic position and its 
characteristic fluctuation. The approximate probability distribution is now 

(2.11) i>((ix, a; k, X) = Z~^ exp |—/3(VbMD — A ‘ ^)| 

Again, the Ai ^re Lagrange multipliers which must satisfy the analog of 

(USD, 

( 2 . 12 ) E^[ai] = a. 


For this case, we can provide expressions for both the partition function and 
the Au The approximate partition function is 

(2.13) i = , Z, = ^exp{-^|x,-XA^jdxi, 

where |D| is the volume of the computational domain. The chemical poten¬ 
tials Au defined by (12.120 . satisfy (see (IA.2I) 1 


(2.14) 


E"[a,] 


D\ + e^i^iZi 


These can then be solved for each i to obtain: 
( 2 . 15 ) = 


The expression for the chemical potentials, (I2.15p . can also be nsed to write 
a simplihed expression for the partition function: 

(2.16) i = 

1 - Cj 

The approximate free energy is now given by 


-/3 ^\ogZ + jl-c 


N 

j3~^ ^ Ci log Ci + {1- Ci) log(l - Ci) 
i=l 


+ {ci - 1) log \D\ - a log Zi. 


(2.17) 
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Note that because D is a bounded set, X* is not exactly the mean 
position of atomic site i. See Appendix El for details. 

2.2.3. Relative Entropy and the Generalized Variational Approach. The ba¬ 
sis for optimizing (I2.10h so that u provides the best match to v is in¬ 
spired by the Gibbs-Bogoliubov inequality and the variational Gaussian ap¬ 
proach [MllIS]. As alluded to in the introduction, this is closely related to 
the relative entropy metric. 

Let us assume that V, the primitive potential in this problem, is bounded 
on the set dB In particular, we assume 

Assumption 2. There exists a constant C > 0 such that for all x G 
and for all a, 

(2.18) |VbMD(x, a)| < C. 

Under these assumptions, we have the necessary absolute continuity to 
proceed with using TZ: 

Proposition 1. Subject to Assumptions E and (H u <C u with Radon- 
Nikodym derivative 

(2.19) ^(x, a) = ^ exp |/3(Udmd - Udmd + /3(A - /^) • a} 
and 

(2.20) 7^(^>||u) =/3E'^[AU]+/3(A-/x)-E^[a] + logZ-logi<oo. 

Proof. Let A be any measurable subset of the state space, ({0,1} x D)^, of 
the form 

( 2 . 21 ) A = Ut^{aixB,) 

where each Bi is a Lebesgue measurable subset of D C If v{A) = 0, 
then 

0 = / exp{-/3(UDMD(x,a) - p • a)}dx > 

where \Bi\ is the Lebesgue measure of the set Bi and C is the constant from 
Assumption [2j Since the pi are assumed to be finite, we can thus infer that 
the only way u(A) = 0 is if at least one of the Bi has Lebesgue measure 
zero, regardless of the particular site occupancy values, a*. Then, since the 
ki defining Udmd are finite, 

z)(A) = Z~^ / exp I -/?(VbMD(x, a) - p • a) [ dx 

( 2 . 22 ) JbixB2...xBn ^ ^ 

< Z-^ef^^i\f^*^Ui\Bi\ = 0 . 

^While this formally excludes potentials which diverge near the origin, such as the 
Lennard-Jones potential, such divergences are routinely avoided by choosing a cut-off at 
an appropriately small interatomic separation. 
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Since this holds for all sets of type (|2.2ip . we conclude it that for all mea¬ 
surable subsets of ({0,1} X D)^ for which iy{A) = 0 we must have u(^) = 0. 
This gives us absolute continuity of the measures. 

The Radon-Nikodym derivative is then given by (I2.19P and TZ is given by 
(12.201) . By our assumptions, TZ will be finite. □ 


Since 7^ > 0, we can use (I2.12|) and (12.171) to express (I2.20p as 

N 

F < E''[AV] -h (3~^'^Ci\ogCi + (1 - Cj)log(l - Cj) - ct logZj 


(2.23) 


i=l 


+ {Ci - l)log|T>|. 


We thus define the approximate DMD free energy, which is an upper bound 
on the true DMD free energy, as 


(2.24) 


N 

F = E'"[AV] -h ^CjlogCi + (1 - Ci)log(l - q) 

i=l 


+ (Ci - l)log|D 


a log Zi 


The DMD algorithm for the vacancy problem proceeds in two steps: 

(1) Find a minimizer of F over k and X, with Xj € D and ki > 0. 

(2) Approximate the dynamics of (j2.8p . substituting F for F. 

Thus, since the Gibbs-Bogliubov inequality is equivalent to the statement 
that the relative entropy is non-negative, the first step in the above algorithm 
is to find the best approximation, with respect to relative entropy, of u over 
a class of distributions of type u. 


2.3. Free Energy and Relative Entropy for a Binary Mixture. Mir¬ 
roring our examination of the vacancy problem, we consider the analogous 
formulation for a binary alloy. 


2.3.1. Canonical Ensemble. For a binary mixture, many of the calculations 
are similar, or even simpler, than for the case of an elemental material with 
vacancies. First, we formulate the uab distribution: 

(2.25) UAB(dx, a) = exp {-/3(Vab - p • a)} dx. 

We continue to assume Assumption [1] holds, adapted to the binary mixture 
case. 


2.3.2. Approximate Ensemble. Now, instead of the approximate potential 
given by (I2.10p for the vacancy problem, we assume 

(2.26) Vab = |xi-X,p. 

i 

The distinction here is that because each site always contains an atom of 
species A or B, the potential is always non-zero. In contrast, the potential 
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associated with a vacancy is always zero and a vacancy is never subject to 
a force. Proceeding with the potential (12.261) . we find 

(2.27) z>AB(dx,a;k,X) = exp |-/3(1 /ab - A-a)|dx, 

where the chemical potentials are chosen to satisfy = Cj. In this 

case 

(2.28) ZAB = n)Ii(l + e^^‘)z„ 

where Zi dehned as in (I2.13p . Hence, we can immediately write the chemical 
potential as 

(2.29) jli = log • 

The free energy is then given by 

N 

(2.30) - \ogZ + n-c = ^CjlogCi + (1 - Ci)log(l - q) - logZ*. 

i=l 


2.3.3. Relative Entropy for the Binary Mixture. We next consider the rela¬ 
tive entropy minimization problem for the binary mixture, requiring uab "C 
UAB- Using the same approach as in the vacancy case, we will also assume 
the boundedness of Vab j in the same spirit as Assumption [2j 


Proposition 2. Under Assumptions {1\ andl^for the binary mixture ease, 
h'AB ^ i^AB with the associated Radon-Nikodym derivative, and 

(2.31) 7e(i>AB||uAB) = /3E''AB[Al/]-h/3(/i-/i)-c-logZAB + log2'AB < oo. 

As before, we can reformulate this as a free energy statement, 

-/3“Mog Zab + /^ • c 

-V-' 

^AB 

< E'^ab [ap] ^ a log Ci + {1- Ci) log(l - Ci) - log Zj. 

i=l 

' -V-' 

Tab 


The simulation now proceeds as above, with a minimizer TZ over k and X 
while the Cj evolve. 


3. Discussion 


3.1. Relation to Existing DMD Formulations. 
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3.1.1. Vacancy Problem. We compare (|2.24p to the DMD free energy for the 
vacancy problem found in [T5l[23l[23]. Indeed, if one rewrites equation (3) 
from the original 2011 DMD paper m for the case of the pair potential 
(in the classical approximation - i.e., for Planck’s constant equals zero), the 
expression is: 


F2011 = E {M-) ^ [M]) ^ fl,, ^(1-^ - -d)e 








|x.-XM 


dxidxj 


( 3 . 1 ) 


iV ^ 

+ /3“^ X) 2 “* ^ + (1 - Ci) log(l - Ci). 


Some of the differences between (13.ip and (I2.24p can be reconciled by the 
use of a mean field approximation, by which (12.21) is replaced with 


(3.2) PDMD,mf — ^ 

i<j 

where the Cj take continuous values between 0 and 1. This approximation 
simplifies some of the computations. For instance, the Hi are now explicit: 

(3.3) Hi = log > 

and the numerical estimation of E'^[IdDMD] is simplified. Notice that (13.3p is 
the same quantity as we obtained in our examination of the binary mixture, 
(12.291) . This will hold generically when the potential for our distribution, 
whether true or approximate, does not explicitly depend on a. Furthermore, 
in the mean field case, the finiteness of the partition function, assumed in 
AssumptionPll is implied by boundedness of the mean field potential over the 
set D for all c G [0,1]'^; thus, an assumption like Assumption [2] is required. 

However, the mean field approximation does not fully account for the 
differences. Part of the discrepancy may be attributed to the choice of the 
state space; Sarkar et al. [I5l[^ use ({0,1} x M'’*)'^. Formally, as D ^ M'’*, 

(3.4) ^/3 ^ ^ Cj , 
and (I2.24P tends to 

N 

(3.5) F^T2oii+/ 3-^ J^(ci-l)log|D|, 

i=l 

where we have not taken \D\ to the limit in the last expression. This de¬ 
pendence on \D\ is a finite size correction to the free energy. As D —>■ 
this becomes unbounded. However, this term does not alter the algorithm 
since we are only concerned with differences of free energies rather than 
absolute magnitudes. Indeed, during the first part of the algorithm, where 
a local minimizer of F is sought over k and X with fixed c, there are no 
contributions oc log \D\. And in the second part of the algorithm, under the 
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dynamics of type ()2.8p . where F is used in place of F, the log |D| terms 
cancel one another. 

The log \ D\ terms also cancel under the more sophisticated dynamics used 
in [U llSlIM] , There, the dynamics are given by the master equation 

(3.6) Cj = |cj(l — — Ci(l — , 

j€N{i) 


with 

(3.7) 


/i = log 

UCi 


1 - a 


Because of the differentiation with respect to c, in the previous expression, 
the log \D\ term does not appear in the differences, fi — fj. 

A more substantive difference between the existing DMD literature and 
our analysis is the notion of the ensemble averaged site positions. As derived 
in the Appendix (|A.3p . 


E" xd = 


(3.8) 


fjj y.idy.i + eddi x* exp |-^ |xj - Xj^j dx* 
\D\ + eddi J^exp |-^ |xj - Xj|^| dx* 

Xidxi Id 1^* — Xj| I 


dx,: 


= (1 - Ci) 


\D\ 


+ Ci 


Z^ 


where we have used (I2.15p to simplify the expression. For simulations in 
computational domains that are symmetric about the origin, such as D = 
(—L, L)'^ or a hypersphere, the first term in this expression vanishes and 




Ci^^i . 


An alternative notion of mean position could be useful in this case. Con¬ 
sider: 


(3.9) 


E^[aiXi] 

E'^lad 


1 


I 2 


= —/ xjexp<j|xi-Xi| 
^iJD 


dxi 


where we have made use of (|A.4p . Now, as D tends to (13.9p recovers 
Xj. This weighted averaging is inspired by the mathematical theory of 
multiphase flow (see, for instance, M)- 

Working with is inherently problematic, as the measures defined as 
in (j2.4p and ()2.1ip do not lead to well defined probability measures, even 
under the mean held approximation. Consider, for example, u with N = 1. 
For this problem the partition function would be 
1 


(3.10) 


E 

ai=0 ‘ 


exp \ -P 


aiki 


x^ X^ I 


+ fdjliai > dxi = oo. 


since, in the case oi = 0, we are integrating dxi over the whole space. 

However, we contend that while (13.111 may give rise to physically con¬ 
sistent simulations, it is not based on a variational principle, and instead, 
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practitioners should use (I2.1ip . In order to make use of TZ and arrive at a 
variational formulation, it is essential that the distributions be well defined 
probability distributions, with finite partition functions and that absolute 
continuity holds. We note that this is not an artifact of our mathemati¬ 
cal analysis based on relative entropy. As the Gibbs-Bogliubov inequality 
is a restatement of the non-negativity of relative entropy, it has the same 
underlying assumptions of absolute continuity and well defined measures. 

One way to correct (j3.8p is to alter the choice of (I2.in[) . Suppose we 
mimic what is done in the binary mixture (which does not suffer from these 
problems), and took 

(3.11) l^DMD = y . 

i 

We would replicate (j2.28j) for the partition function and, with this revised 
value, the free energy would be 


N 


2=1 


(3.12) Fdmd = /3 ^ log Cj(1 - Ci) log(l - Cj) - log Zj. 

Now, as D ^ M'^, we obtain 

.AdMD ~ ^ ^ Cj log Ci -I- (1 — Ci) log(l — Ci) 


N 


(3.13) 


2 = 1 


2 I ^ 27r 


This formulation has the advantage that E^[xj] ^ Xj as H ^ However, 
notice now that there is no Cj multiplying the last expression in (j3.13p . as 
in (j3.ip . 


3.1.2. Binary Mixture Problem. As noted, if we make a mean field approx¬ 
imation, then as D —)> we recover the expression found in [3], which we 
do not reproduce here. 


3.2. Advantages of Relative Entropy. One of the main advantages of 
the formulation given here is that the question of relative entropy mini¬ 
mization is a rigorously defined variational problem, forcing us to confront 
problems such as that associated with domain size. Indeed, for DMD we 
have 

Theorem 1. Given K G (0, oo) and the open bounded subset D let 

(3.14) Ad,k = {7 of the form (I2.11h | Xj € D, ki & [0; AT]} . 

Then ifi'jn} is a minimizing sequence ofTZ{-\\v), it has a subsequential limit 
such that 
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and 

(3.16) 


TV 

7n ^ 7*- 

Note that by our previous assumptions on V and /x, that they are bounded 
on D for Cj G (0,1), we are assured that the elements of Ad,k are absolutely 
continuous with respect to such a u, and have finite Dki- 

Proof. We claim Ad,k is weakly compact; see Appendix [Bj Therefore, any 
sequence has a weak limit, 7 ^^, ^ 7*. Then, by the lower semicontinuity 
of 7^(-||u) (Proposition 2.1 of [20]) we have (13.151) . then we can infer (I3.16P 
(Lemma 2.4 of [2Q|). □ 

Thus, the free energy minimization part of the algorithm is well posed in 
the sense that if we take a minimizing sequence it has a subsequential limit. 

Another advantage of this formulation of the problem is that it allows one 
to consider more general parameterizations of the synthetic distribution, u. 
Indeed, one could imagine any distribution parameterized by some collection 
of variables, denoted collectively by p, and proceed as above; first, one 
finds a local minimizer of 7^(z)p||u). Using that, the free energy gradients 
with respect to c are computed at this stationary point, and this drives the 
dynamic evolution of the Cj’s. Provided this admissible class is closed (in 
some sense), we are ensured that its minimization problem is well posed too. 


3.3. Open Problems and Future Work. Here, we have presented a 
mathematical framework for DMD, and we have given particular attention 
to the question of free energy minimization. One point we have not ad¬ 
dressed is the temporal evolution problem, and how the master equation 
will be influenced by the use of the approximate free energy in place of the 
true DMD free energy. Indeed, we note that while the minimization of free 
energy ensures that u is close to u in the sense of relative entropy (and hence, 
total variation), we cannot conclude that 


(3.17) 


dF dP 

dci dci ’ 


errors in the dynamics of both (12.Sp and (13.6p are small. As of now, this 
remains an unconstrained approximation that merits investigation. It can 
be shown that 

dF 

(3.18) — = 

and, at a minimizer of IZ, 


(3.19) 


dP 

dci 


fii + (5 Covp (AU, dcifi ■ a). 


Calculations of these can be found in the appendix. Similar expressions hold 
under the mean field approximation. 
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The validity of the mean field approximation used in [Illl5tl2^ remains 
to be explored. As noted, this dramatically simplifies the calculations and 
immediately gives the “true” chemical potentials (13.3p . The solvability in 
the general case is also an open problem. 

Finally, there is the question of how, or in what sense, does DMD, as 
a model, approximate any specific, more primitive MD model, such as 
Langevin dynamics. More specifically, since DMD removes the stochas¬ 
tic nature of primitive MD models, replacing the displacive dynamics with 
a finite-temperature, energy minimization, DMD is deterministic. Simi¬ 
larly, the stochastic or random walk-like diffusive dynamics in primitive MD 
becomes deterministic in DMD. The question of when this is and is not 
acceptable remains to be explored. 


Appendix A. Detailed Calculations 


A.l. Partition Functions. To obtain the approximate partition function 
(j2.13j) . we employ the following procedure: 


Z = ^ J exp |—/?(l/(x, a; k, X) —/i • a 


dx 


(A.l) 


— '-H=\ 


exp |xi - Xil^-h/3/iiai| dxj 

|-^ |xj - Xjpj dxjj 


= |D| + e^^/ exp -^|x,-X 


= 


,[\D\ + e^^Z,]. 


A.2. Chemical Potentials. To obtain (I2.15p . we apply (lA.ip to obtain 
(A.2) 


lE'K] = 




Z 




1,^ -V" I 

2 I I 




dx,; 




\D\ + e^i^^Zi 

Since E^[aj] = c,, we solve for fli in terms of Cj 


A.3. Mean Position. The expectation value of the atomic position x* is 
(A.3) 

Z 


xd = 




E 


Xj exp 


D 


1 


|D| -|- ed/'iZ, 


V 

Xjdxj -h j Xjexp |-^ |xj - Xj|^| dxj^ 


r ^ajkj I -y- |2 I 

\ 2 1^* + 




dxi 


i \J D 
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and the weighted mean position is given by 


(A.4) 


E‘'[aiXi] 


Z 



I Xi 



dxi 


1 


\D\ + 


,/iii j X • exp I 



A.4. Mean Potential. The expectation value of the mean potential I^md 
is obtained as follows: 


(A.5) 


N 


E"[1>dmd] = 

i=i 


i„ V |2 

./V T .Z\. T I 

2 I 1 


Z 


i=i 

i=i 


e/3Ai 


Y l^i - Xjf exp I \xj - Xj |2| dxj 


2Zj jd 


|xj - Xjl^ exp |-^ |xj - Xjf} dxj 


A.5. Free Energy Gradients. The DMD free energy gradient, with re¬ 
spect to Ci, is computed as follows 


dF 

dci 


I 1 dZ dfi 


-/3- 




Qfj/ 

-/3Vdmd + fJ- ■ a} dx'f + — ■ c + Hi 




c -h /li = Hi- 


This is (j3.18p . The gradient of F with respect to any of the parameters 
defining u vanishes at the minimizer of F. To see this, recall that since we 
have a minimizer of IZ, 


c>7^(^>||u) 

dp 


= 0 


for any parameter, p, such as ki and Xj. Since we can write 


l3-^n = F-F 
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and F does not depend on the parameters, when a derivative is taken with 
respect to Cj, chain rule terms involving the parameters vanish at the mini- 
mizers. Therefore, 

^ - A-ilM M '■ 

dci Z dci ^ dci Qci 

= /Ij + i ^ f AV(3{dcip- ■ exp i-(3V + /?/i • a| 

- ^¥F[AV] J /3{dcifi ■ a) exp |-/3y + /3A • a}^ , 
and this gives us (|3.19|) . 


Appendix B. Weak Compactness of the set of Measures 


Lemma 1. The set Ad,k is weakly compact, in the sense that if {7n} is 
any sequence in Ad,k, it has a weakly converging subsequence in Ad,k- 

Proof. Given any sequence of 7 „ € Ad,k, we have seqnences ^ 

X [0, Since the set is compact, it has a convergent subsequence, 

xG-) ^ ^ kl*'> 

Let 7 * be the measure associated with and k^*). Clearly, 7 * G Ad,k- 
Given any bounded continuons function / on the set X, we will now show 

From (I2.16p . so long as the Cj € (0,1) = Z*. For brevity, let 

VbMD(x,a;X("-),k("-)) = l/(”^)(x,a) ^ yW(x,a) 

Therefore, when we take differences, 

< 

< 

Since x and a are elements of bounded sets, and 

g-/3(U(™) (x,a)-UW (x,a))+/3(/xM -/i(*))-a 


/ l/(x,a)| 


Id^ 


g-/3U(™)(x,a)+,3A^’"'-a _ g-^UW(x,a)+/3A^*^-f 


Y / 


(x,a)+/3A^*^-i3 




depends continnonsly upon and \A^\ we then have there exists some 
constant such that for all x € A and a G {0,1}, 


^_^(FM(x,a)-U(*)(x,a))+/3(/x("‘)-/xW)-a _ 


< c 



xW + 
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Therefore, 

IET""- [/] - E'^*[/]| < C - xW + - kW ) E>[|/|] 

and we have weak convergence. 


□ 
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